library(CellChat)
library(patchwork)
library(tidyverse)
library(Seurat)
library(tidyseurat)
library(ggpubr)

# initial CellChat object --------
# normalized data matrix is needed
crc4 <- read_rds('CRC-I/data/crc_merge4_immune.rds')

liu23 <- crc4 |>
  filter(dataset == 'liu2023')

# a dataframe with rownames containing cell mata data
meta <- liu23@meta.data

# directly from Seurat or SCE
obcc <- createCellChat(liu23, group.by = 'manual_main')

obcc |> glimpse()

## set ppi db ---------
showDatabaseCategory(CellChatDB.human)

glimpse(CellChatDB.human)
CellChatDB.use <- subsetDB(CellChatDB.human, search = "Secreted Signaling")

obcc@DB <- CellChatDB.use

# identify genes in db
obcc <- subsetData(obcc)

# if use multisession, no progress bar is shown
# and improve in speed is not significant
# future::plan("multisession", workers = 4)

obcc <- identifyOverExpressedGenes(obcc)
obcc <- identifyOverExpressedInteractions(obcc)

## infer cell-cell network ------------
obcc <- computeCommunProb(obcc) # cost ~3m

## use looser to find weaker expressed gene
## default triMean method equals ~.25 trim
obcc1 <- computeCommunProb(obcc, type = 'truncatedMean', trim = .1)

# try signaling of interest to find proper threshold
obcc |> computeAveExpr(features = 'FCGR2B', type = 'truncatedMean', trim = .05)

## or to consider cell prop
## better for no-sort-enriched sc data
obcc <- computeCommunProb(obcc, population.size = TRUE)

## filter out interactions with low cell number in one group
obcc <- filterCommunication(obcc, min.cells = 10)

# extract result as a data.frame
## default ligand-pair level
res_obcc <- subsetCommunication(obcc)

## or pathway level
res_obcc_pathw <- subsetCommunication(obcc, slot.name = 'netP')

res_obcc |> head()
res_obcc_pathw |> head()

# specify source & target in interaction
subsetCommunication(obcc, sources.use = c(1,2), targets.use = c(4,5))

# specify signaling module
subsetCommunication(obcc, signaling = 'TGFb')

# aggregate network and visualize them
obcc <- aggregateNet(obcc)

groupSize <- table(obcc@idents) |> as.numeric()

# graph param: Map From ROW: 1 row * 2 col, expanded
par(mfrow = c(1,2), xpd=TRUE)
obcc@net$count |>
  netVisual_circle(vertex.weight = groupSize,
                   weight.scale = T,
                   label.edge= F,
                   title.name = "Number of interactions")

obcc@net$weight |>
  netVisual_circle(vertex.weight = groupSize,
                   weight.scale = T,
                   label.edge= F,
                   title.name = "Interaction weights/strength")

## split view from each cell type
mat <- obcc@net$weight

par(mfrow = c(2,4), xpd=TRUE)
for (i in 1:nrow(mat)) {
  mat2 <- matrix(0, nrow = nrow(mat), ncol = ncol(mat), dimnames = dimnames(mat))
  mat2[i, ] <- mat[1, ]
  mat2 |>
    netVisual_circle(vertex.weight = groupSize,
                     weight.scale = T,
                     edge.weight.max = max(mat),
                     title.name = rownames(mat)[i])
}

## visualize network ----------
pathways.show <- c("TNF") 

## Hierarchy plot
# define `vertex.receiver` to specify clusters in left portion of plot 
obcc |>
  netVisual_aggregate(signaling = pathways.show,
                      vertex.receiver = c(1:2,4),
                      layout = 'hierarchy')
## Circle plot
par(mfrow=c(1,1))

obcc |>
  netVisual_aggregate(signaling = pathways.show, layout = "circle")

## Chord diagram
par(mfrow=c(1,1))
obcc |>
  netVisual_aggregate(signaling = pathways.show, layout = "chord")

## Heatmap
par(mfrow=c(1,1))
obcc |>
  netVisual_heatmap(signaling = pathways.show, color.heatmap = 'Reds')

## fine-tune Chord diagram
levels(obcc@idents)

group.cellType <- c('S','S','I','S','I','I','I') |>
  set_names(levels(obcc@idents))

obcc |>
  netVisual_chord_cell(signaling = pathways.show,
                       group = group.cellType,
                       title.name = paste0(pathways.show, " signaling network"))

## show pairs in pathway
pairLR.CCL <- obcc |>
  netAnalysis_contribution(signaling = pathways.show, return.data = T)

pairLR.CCL$gg.obj

## show most enriched pair
LR.show <- pairLR.CCL$LR.contribution[1,] |> rownames() # show one ligand-receptor pair
# Hierarchy plot
obcc |>
  netVisual_individual(signaling = pathways.show,
                       pairLR.use = LR.show,
                       vertex.receiver = c(1,2,4),
                       layout = 'hierarchy')
# Circle plot
obcc |>
  netVisual_individual(signaling = pathways.show, pairLR.use = LR.show)

obcc |>
  netVisual_individual(signaling = pathways.show, pairLR.use = LR.show, layout = 'chord')

# save all plot in one time!
obcc <- computeCommunProbPathway(obcc)
pathways.show.all <- obcc@netP$pathways
# check the order of cell identity to set suitable vertex.receiver
levels(obcc@idents)
vertex.receiver = c(1,2,4)
for (i in 1:length(pathways.show.all)) {
  # Visualize communication network associated with both signaling pathway and individual L-R pairs
  obcc |> netVisual(signaling = pathways.show.all[i],
                    vertex.receiver = vertex.receiver,
                    layout = "hierarchy")
  # Compute and visualize the contribution of each ligand-receptor pair to the overall signaling pathway
  gg <- obcc |>
    netAnalysis_contribution(signaling = pathways.show.all[i])
  paste0(pathways.show.all[i], "_L-R_contribution.pdf") |>
    ggsave(plot=gg, width = 3, height = 2, units = 'in')
}

# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
obcc |>
  netVisual_bubble(sources.use = 4,
                   targets.use = c(3,5:7),
                   remove.isolate = FALSE) +
  theme_pubr(x.text.angle = 45,legend = 'right')

# show all the significant interactions (L-R pairs) from some cell groups (defined by 'sources.use') to other cell groups (defined by 'targets.use')
# the width of chord indicate strength of interaction
obcc |>
  netVisual_chord_gene(sources.use = 4,
                       targets.use = c(3,5:7))

# show chords of pathways
obcc |>
  netVisual_chord_gene(sources.use = 4,
                       targets.use = c(3,5:7),
                       slot.name = 'netP')

# show signaling genes in vlnplot
obcc |>
  plotGeneExpression(signaling = 'CCL')

# Compute the network centrality scores ------
obcc <- obcc |>
  netAnalysis_computeCentrality(slot.name = "netP")

# Visualize the computed centrality scores using heatmap, allowing ready identification of major signaling roles of cell groups
# default size can be too small
obcc |>
  netAnalysis_signalingRole_network(signaling = 'CCL')

# scatter plot show incoming & outgoing strength of cell types
obcc |>
  netAnalysis_signalingRole_scatter(signaling = 'CCL')

# show IO pathways in heatmap
ht1 <- obcc |> netAnalysis_signalingRole_heatmap(pattern = "outgoing")
ht2 <- obcc |> netAnalysis_signalingRole_heatmap(pattern = "incoming")
ht1 + ht2

# identify global patterns ----------
library(NMF)
library(ggalluvial)

## outgoing patterns, can take a long time
selectK(obcc, pattern = 'outgoing')
# choose k number when both Cophenetic and Silhouette values begin to drop suddenly
nPatterns = 6
obcc <- obcc |>
  identifyCommunicationPatterns(pattern = "outgoing", k = nPatterns)

# river plot
obcc |> netAnalysis_river(pattern = "outgoing")

# dot plot
obcc |> netAnalysis_dot(pattern = "outgoing")

## incoming patterns
selectK(obcc, pattern = 'incoming')
# optimal k number can be different for I & O
nPatterns = 2
obcc <- obcc |>
  identifyCommunicationPatterns(pattern = "incoming", k = nPatterns)

# river plot
obcc |> netAnalysis_river(pattern = "incoming")

# dot plot
obcc |> netAnalysis_dot(pattern = "incoming")

# identify pathways with similar sender & receiver ------------ 
obcc <- obcc |> computeNetSimilarity(type = "functional")
obcc <- obcc |> netEmbedding(type = "functional", umap.method = 'uwot')

obcc <- obcc |> netClustering(type = "functional", do.parallel = FALSE)

# Visualization in 2D-space
obcc |> netVisual_embedding(type = "functional", label.size = 3.5)

# structural similary only consider network geometry, agnostic of node identity
obcc <- obcc |> computeNetSimilarity(type = "structural")
obcc <- obcc |> netEmbedding(type = "structural", umap.method = 'uwot')

obcc <- obcc |> netClustering(type = "structural", do.parallel = FALSE)

# Visualization in 2D-space
obcc |> netVisual_embedding(type = "structural", label.size = 3.5)

# Compare data for same composition ----------
liu23 <- liu23 |>
  filter(!str_detect(manual_main, 'Dividing'))

l23tt <- liu23 |>
  filter(genotype == 'TT')

l23ctl <- liu23 |>
  filter(genotype != 'TT')

## identify interaction respectively ------------
obcctt <- l23tt |>
  createCellChat(group.by = 'manual_main')

obccctl <- l23ctl |>
  createCellChat(group.by = 'manual_main')

showDatabaseCategory(CellChatDB.human)

CellChatDB.use <- subsetDB(CellChatDB.human, non_protein = FALSE)

obcctt@DB <- CellChatDB.use
obccctl@DB <- CellChatDB.use

# identify genes in db
obcctt <- obcctt |>
  subsetData() |>
  identifyOverExpressedGenes() |>
  identifyOverExpressedInteractions() |>
  computeCommunProb(population.size = TRUE) |>
  aggregateNet() |>
  computeCommunProbPathway() |>
  netAnalysis_computeCentrality()

obccctl <- obccctl |>
  subsetData() |>
  identifyOverExpressedGenes() |>
  identifyOverExpressedInteractions() |>
  computeCommunProb(population.size = TRUE) |>
  aggregateNet() |>
  computeCommunProbPathway() |>
  netAnalysis_computeCentrality()

## merge 2 processed data ---------
obcc <- list(obccctl, obcctt) |>
  mergeCellChat(add.names = c('II+IT', 'TT'))

### compare total interaction number & strength --------
g1 <- obcc |>
  compareInteractions(group = c('II+IT', 'TT'),
                      color.use = c('blue','red'))

g2 <- obcc |>
  compareInteractions(group = c('II+IT', 'TT'),
                      color.use = c('blue','red'), measure = 'weight')

g1 + g2 & theme_pubr(legend = 'none')

### diff among net ---------
### chord plot
par(mfrow = 1:2, xpd = T)
obcc |> netVisual_diffInteraction(weight.scale = T)
obcc |> netVisual_diffInteraction(weight.scale = T, measure = 'weight')

### heatmap
obcc |> netVisual_heatmap()
obcc |> netVisual_heatmap(measure = 'weight',
                          font.size = 12,
                          font.size.title = 12)

### narrower chord
group.cellType <- c(rep('CD4T',4),rep('CD8T',4),rep('Myeloid',4)) |>
  factor()

levels(group.cellType)

obcc_list <- list(obccctl, obcctt) |>
  map(\(x)mergeInteractions(x, c('CD4T','CD8T','Myeloid'))) |>
  set_names(c('II+IT','TT'))

obcc3 <- obcc_list |>
  mergeCellChat(add.names = c('II+IT','TT'))

weight.max <- obcc_list |>
  getMaxWeight(slot.name = c('idents','net','net'),
               attribute = c('idents','count','count.merged'))

#### number of interactions among 3
par(mfrow = 1:2, xpd = T)
for (i in 1:2) {
  obcc_list[[i]]@net$count.merged |>
    netVisual_circle(weight.scale = T, label.edge = T,
                     edge.weight.max = weight.max[3],
                     edge.width.max = 12,
                     title.name = str_c('Number of interactions - ',
                                        names(obcc_list)[i]))
}

#### differential number & strength
par(mfrow = 1:2)
obcc3 |>
  netVisual_diffInteraction(weight.scale = T,
                            measure = 'count.merged', label.edge = T)

obcc3 |>
  netVisual_diffInteraction(weight.scale = T,
                            measure = 'weight.merged', label.edge = T)

## incoming & outgoing --------------
num.link <- obcc_list |>
  sapply(\(x){
    y <- x@net$count
    rowSums(y) + colSums(y) - diag(y)
  })

weight.min.max <- c(min(num.link), max(num.link))  

gg <- list()

for (i in 1:2) {
  gg[[i]] <- obcc_list[[i]] |>
    netAnalysis_signalingRole_scatter(title = names(obcc_list)[i],
                                      weight.MinMax = weight.min.max)
}

wrap_plots(gg)

### identify specific signal in one cell type
obcc |>
  netAnalysis_signalingChanges_scatter(idents = 'Myeloid',font.size = 12,
                                       color.use = c('grey','blue','red')) +
  theme_pubr(legend = 'right')

## singaling network structure ------------
### functional (default) -----
obcc <- obcc |>
  computeNetSimilarityPairwise() 

obcc |>
  netEmbedding(pathway.remove = 1, n_neighbors = 8) 

obcc |>
  netClustering()

Similarity <- methods::slot(obcc, 'netP')$similarity[['functional']]$matrix[['1-2']]

obcc |>
  computeNetSimilarityPairwise() |>
  rankSimilarity()

### information flow of pathways ---------
gg1 <- obcc |>
  rankNet(do.stat = T, color.use = c('blue','red'),font.size = 10) +
  NoLegend()

gg2 <- obcc |>
  rankNet(do.stat = T, color.use = c('blue','red'), font.size = 10, stacked = T)

gg1 + gg2

### outgoing pattern ----------
pathway.union <- union(obcc_list[[1]]@netP$pathways,
                       obcc_list[[2]]@netP$pathways)

ht1 <- obcc_list[[1]] |>
  netAnalysis_signalingRole_heatmap(pattern = "outgoing",
                                    signaling = pathway.union,
                                    title = names(obcc_list)[1],
                                    width = 5, height = 6)

ht2 <- obcc_list[[2]] |>
  netAnalysis_signalingRole_heatmap(pattern = "outgoing",
                                    signaling = pathway.union,
                                    title = names(obcc_list)[2],
                                    width = 5, height = 6)

#ComplexHeatmap::draw(ht1 + ht2, ht_gap = unit(0.5, "cm"))

### differential ligand pairs -----------
obcc |>
  netVisual_bubble(angle.x = 45,
                   sources.use = 5,
                   targets.use = 2:3,
                   comparison = 1:2,
                   color.text = c('blue','red'))

obcc |>
  netVisual_bubble(angle.x = 45,
                   sources.use = 5,
                   targets.use = 2:3,
                   comparison = 1:2,
                   color.text = c('blue','red'),
                   max.dataset = 2,
                   remove.isolate = T,
                   title.name = 'Increased signaling in TT')


obcc |>
  netVisual_bubble(angle.x = 45,
                   sources.use = 5,
                   targets.use = 2:3,
                   comparison = 1:2,
                   color.text = c('blue','red'),
                   max.dataset = 1,
                   remove.isolate = T,font.size = 12,
                   title.name = 'Decreased signaling in TT')

#### use DEA method ---------
obcc <- obcc |>
  identifyOverExpressedGenes(group.dataset = "datasets",
                             pos.dataset = 'TT',
                             features.name = 'TT.merged',
                             only.pos = FALSE, thresh.pc = 0.1,
                             thresh.fc = 0.05,thresh.p = 0.05,
                             group.DE.combined = T)

net <- obcc |>
  netMappingDEG(features.name = 'TT.merged',
                variable.all = T)

net.up <- obcc |>
  subsetCommunication(net = net, datasets = "TT",
                      ligand.logFC = 0.05, receptor.logFC = NULL) |>
  extractGeneSubsetFromPair(obcc)

net.down <- obcc |>
  subsetCommunication(net = net, datasets = "II+IT",
                      ligand.logFC = -0.05, receptor.logFC = NULL) |>
  extractGeneSubsetFromPair(obcc)

